   // StepWaveProb.C

#include <iostream>
#include <fstream>
#include <string.h>
#include <math.h>

#include "NOX_Common.H"
#include "Epetra_SerialComm.h"
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_Import.h"
#include "Epetra_CrsGraph.h"
#include "Epetra_CrsMatrix.h"

#include "Global.H"
#include "StepWaveProb.H"


using namespace std;


StepWaveProb::StepWaveProb(int m, int n, double p0, double p1, double p2, double a1, double a2, double a3)
{
  M = m;
  N = n;
  P0 = p0;
  P1 = p1;
  P2 = p2;
  A1 = a1;
  A2 = a2;
  A3 = a3;
  firstFill = true;

  Epetra_SerialComm comm;
  map = new Epetra_Map(M*N, 0, comm);
  initial = new Epetra_Vector(*map);
  J = new Epetra_CrsMatrix(Copy, *map, 9);

  double lambda = computeLambda();
  initQ(lambda);
  computeInitial(lambda);
}

StepWaveProb::StepWaveProb(char *solfile)
{
  firstFill = true;

  ifstream fin(solfile); if (fin.fail()) { cerr << "Failed to open initial solution file, exiting" << endl; exit(1); }
  char header[32];
  fin >> header; if (strcmp(header, "step")) { cerr << "Error: initial solution file provided does not begin with \"step\", exiting" << endl; exit(1); }
  fin >> M;
  fin >> N;
  fin >> P0;
  fin >> Q;
  fin >> P1;
  fin >> P2;
  fin >> A1;
  fin >> A2;
  fin >> A3;

  Epetra_SerialComm comm;
  map = new Epetra_Map(M*N, 0, comm);
  initial = new Epetra_Vector(*map);
  J = new Epetra_CrsMatrix(Copy, *map, 9);

  double v;
  for (int i=0; fin >> v; i++)
    (*initial)[i] = v;
  fin.close();
}

StepWaveProb::~StepWaveProb()
{
  delete map;
  delete initial;
  delete J;
}

double StepWaveProb::computeLambda()
{
  	double gammaMin = 0.0;
  	double hp = fabs(P0) / M;
  	double intGamma = 0.0;
	
	for (int r=0; r<M; ++r)
	{
    	intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp/2;
    	if (intGamma < gammaMin) { gammaMin = intGamma; }
    }  
    
	return newton(-2*gammaMin+1e-2);
}

void StepWaveProb::initQ(double lambda)
{
  double hp = fabs(P0) / M;
 double intGamma = 0.0;
  double next;
  Q = 0.0;
  cerr << "  lambda = " << lambda << endl;
  for (int r=0; r<M-1; ++r)
  {
    intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp/2;
    Q -= (hp/2)/sqrt(lambda + 2*intGamma);
    next = intGamma - (Gamma((r+1)*hp) + Gamma((r+2)*hp))*hp/2;
    Q -= (hp/2)/sqrt(lambda + 2*next);
  }
  int r = M-1;
  intGamma -= (Gamma(r*hp)+Gamma((r+1)*hp))*hp/2;
  Q -= hp/sqrt(lambda + 2*intGamma);
  		
  
  Q = -2*g*Q + lambda;
  //Q = 21.7115;
  cerr << endl << "Starting search for Q* at Q = " << Q << "..." << endl; 
}

void StepWaveProb::computeInitial(double lambda)
{
  double hp = fabs(P0)/M;
  double t2 = (Q - lambda)/(2*g);
  double integral = 0;
  double intGamma = 0;
  double next;
  for (int r=0; r<M-1; ++r)
  {
    intGamma -= (Gamma(r*hp)+ Gamma((r+1)*hp))*hp/2;
    integral -= (hp/2)/sqrt(lambda + 2*intGamma);
    next = intGamma - (Gamma((r+1)*hp)+Gamma((r+2)*hp))*hp/2;
    integral -= (hp/2)/sqrt(lambda + 2*next);
    for (int c=0; c<N; ++c)
      (*initial)[r*N + c] = integral + t2;
  }
  int r = M-1;
  intGamma -= (Gamma(r*hp)+ Gamma((r+1)*hp))*hp/2;
  integral -= hp/sqrt(lambda + 2*intGamma);
  for (int c=0; c<N; ++c)
    (*initial)[r*N + c] = integral + t2;
}
